#include "fft_radix7.h"

#include "fft_common_use.h"

vect_complex_d fft_radix_7(vect_complex_d& f)
{
    int N = f.size();
    vect_complex_d F(N);
    if (N == 7)
    {
        F[0] = f[0] + f[1] + f[2] + f[3] + f[4] + f[5] + f[6];
        F[1] = f[0] + complex_d(0.6234898018587336, -0.7818314824680298) * f[1] + complex_d(-0.22252093395631434, -0.9749279121818236) * f[2] + complex_d(-0.900968867902419, -0.43388373911755823) * f[3] + complex_d(-0.9009688679024191, 0.433883739117558) * f[4] + complex_d(-0.2225209339563146, 0.9749279121818236) * f[5] + complex_d(0.6234898018587334, 0.7818314824680299) * f[6];
        F[2] = f[0] + complex_d(-0.22252093395631434, -0.9749279121818236) * f[1] + complex_d(-0.9009688679024191, 0.433883739117558) * f[2] + complex_d(0.6234898018587334, 0.7818314824680299) * f[3] + complex_d(0.6234898018587336, -0.7818314824680298) * f[4] + complex_d(-0.900968867902419, -0.43388373911755823) * f[5] + complex_d(-0.2225209339563146, 0.9749279121818236) * f[6];
        F[3] = f[0] + complex_d(-0.900968867902419, -0.43388373911755823) * f[1] + complex_d(0.6234898018587334, 0.7818314824680299) * f[2] + complex_d(-0.22252093395631434, -0.9749279121818236) * f[3] + complex_d(-0.2225209339563146, 0.9749279121818236) * f[4] + complex_d(0.6234898018587336, -0.7818314824680298) * f[5] + complex_d(-0.9009688679024191, 0.433883739117558) * f[6];
        F[4] = f[0] + complex_d(-0.9009688679024191, 0.433883739117558) * f[1] + complex_d(0.6234898018587336, -0.7818314824680298) * f[2] + complex_d(-0.2225209339563146, 0.9749279121818236) * f[3] + complex_d(-0.22252093395631434, -0.9749279121818236) * f[4] + complex_d(0.6234898018587334, 0.7818314824680299) * f[5] + complex_d(-0.900968867902419, -0.43388373911755823) * f[6];
        F[5] = f[0] + complex_d(-0.2225209339563146, 0.9749279121818236) * f[1] + complex_d(-0.900968867902419, -0.43388373911755823) * f[2] + complex_d(0.6234898018587336, -0.7818314824680298) * f[3] + complex_d(0.6234898018587334, 0.7818314824680299) * f[4] + complex_d(-0.9009688679024191, 0.433883739117558) * f[5] + complex_d(-0.22252093395631434, -0.9749279121818236) * f[6];
        F[6] = f[0] + complex_d(0.6234898018587334, 0.7818314824680299) * f[1] + complex_d(-0.2225209339563146, 0.9749279121818236) * f[2] + complex_d(-0.9009688679024191, 0.433883739117558) * f[3] + complex_d(-0.900968867902419, -0.43388373911755823) * f[4] + complex_d(-0.22252093395631434, -0.9749279121818236) * f[5] + complex_d(0.6234898018587336, -0.7818314824680298) * f[6];
        return F;
    }
    int limit = N / 7;
    vect_complex_d f_1(limit), f_2(limit), f_3(limit),
        f_4(limit), f_5(limit), f_6(limit), f_7(limit);
    for (int r = 0; r < limit; r++)
    {
        f_1[r] = f[7 * r];
        f_2[r] = f[7 * r + 1];
        f_3[r] = f[7 * r + 2];
        f_4[r] = f[7 * r + 3];
        f_5[r] = f[7 * r + 4];
        f_6[r] = f[7 * r + 5];
        f_7[r] = f[7 * r + 6];
    }
    vect_complex_d F_1, F_2, F_3,
        F_4, F_5, F_6, F_7;
    F_1 = fft_radix_7(f_1);
    F_2 = fft_radix_7(f_2);
    F_3 = fft_radix_7(f_3);
    F_4 = fft_radix_7(f_4);
    F_5 = fft_radix_7(f_5);
    F_6 = fft_radix_7(f_6);
    F_7 = fft_radix_7(f_7);
    for (int k = 0; k < limit; k++)
    {
        complex_d W1(calculate_W(N, k)),
            W2(calculate_W(N, 2 * k)),
            W3(calculate_W(N, 3 * k)),
            W4(calculate_W(N, 4 * k)),
            W5(calculate_W(N, 5 * k)),
            W6(calculate_W(N, 6 * k));
        F[k] = F_1[k] + W1 * F_2[k] + W2 * F_3[k] + W3 * F_4[k] + W4 * F_5[k] + W5 * F_6[k] + W6 * F_7[k];
        F[k + limit] = F_1[k] + w71 * W1 * F_2[k] + w72 * W2 * F_3[k] + w73 * W3 * F_4[k] + w74 * W4 * F_5[k] + w75 * W5 * F_6[k] + w76 * W6 * F_7[k];
        F[k + 2 * limit] = F_1[k] + w72 * W1 * F_2[k] + w74 * W2 * F_3[k] + w76 * W3 * F_4[k] + w71 * W4 * F_5[k] + w73 * W5 * F_6[k] + w75 * W6 * F_7[k];
        F[k + 3 * limit] = F_1[k] + w73 * W1 * F_2[k] + w76 * W2 * F_3[k] + w72 * W3 * F_4[k] + w75 * W4 * F_5[k] + w71 * W5 * F_6[k] + w74 * W6 * F_7[k];
        F[k + 4 * limit] = F_1[k] + w74 * W1 * F_2[k] + w71 * W2 * F_3[k] + w75 * W3 * F_4[k] + w72 * W4 * F_5[k] + w76 * W5 * F_6[k] + w73 * W6 * F_7[k];
        F[k + 5 * limit] = F_1[k] + w75 * W1 * F_2[k] + w73 * W2 * F_3[k] + w71 * W3 * F_4[k] + w76 * W4 * F_5[k] + w74 * W5 * F_6[k] + w72 * W6 * F_7[k];
        F[k + 6 * limit] = F_1[k] + w76 * W1 * F_2[k] + w75 * W2 * F_3[k] + w74 * W3 * F_4[k] + w73 * W4 * F_5[k] + w72 * W5 * F_6[k] + w71 * W6 * F_7[k];
    }
    return F;
}
